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SUMMARY 


The purpose of this paper is to develop a singularity-free analytical solu- 
tion to the J 2 perturbation problem in the Earth satellite theory. Elements 
similar to the classical Poincare elements, but which are in the extended phase 
space, are used. The solution, which has the long-period terms eliminated, is 
developed into a computer algorithm which is then compared for several typical 
examples with numerically integrated results. The error of the analytical solu- 
tion is of the order of 500 to 1000 meters after 100 periods of revolution. 


INTRODUCTION 


The objective of this paper is to develop a singularity-free (for vanishing 
eccentricity and inclination) solution of the J 2 problem in the satellite the- 
ory. The procedure resembles the one described by Lyddane (ref. 1), which 
rederives Brouwer's satellite theory (ref. 2), by using canonical Poincare^ ele- 
ments instead of the Delaunay elements. A comparable procedure is performed in 
this report for the newly developed satellite theory, which uses elements simi- 
lar to the Delaunay elements (DS^) presented by Scheifele in reference 3 and the 
Poincar^ elements (PS<j>) presented by Scheifele and Graf in reference 4. The ca- 
nonical von Zeipel approach is used to arrive at a perturbation solution, which 
will include only the short-period effects due to J 2 . 

1. The fast variable (here canonical true anomaly) is eliminated by the 
use of a canonical transformation, leading to a set of "mean" elements. 

2. The remaining Hamiltonian system, which has a Hamiltonian that is free 
of angle variables, is integrated by quadrature. 

This procedure, when used with the DS(j) and PS<t> elements, necessitates 
changes in both (1) and (2) as follows: 

1 . The generating function that leads to the "mean" elements and is known 
from the DS(1) perturbation theory is reformulated in terms of PS<t> variables. 
Because of the nature of these variables, there will be no more singularities 
when taking the partial derivatives of this new generating function with respect 
to the PS<1> variables. 

2. The solution of the first-order Hamiltonian system in the DS<|> vari- 
ables is rearranged in such a way that it can be expressed entirely in PS<|) 
variables. 


The following observations can be made: 

1 . The use of many convenient regular and carefully selected abbreviations 
will considerably reduce computation time and core requirements. 

2. Typical accuracy is 200 to 500 meters after 50 revolutions when 
compared with numerically integrated solutions. 

3. Although the theory is not restricted to short-period J2 effects 

only, for practical purposes it seems adequate at present to take only the 
short-period J2 solution into account because the higher zonal harmonics 
short-period terms are of the order only. Also, the long-period effects 

due to J2 are the same order as the higher zonal harmonics. 


The author gratefully acknowledges the assistance of Gerhard Scheifele of 
ACM, Inc. for suggesting that this problem be solved and for patiently providing 
invaluable aid in the complicated disciplines of Hamilton mechanics and von 
Zeipel formalism. Dr. Scheifele is the inventor of the Delaunay-Similar and 
Poincar4-Similar elements. He holds credit for first solving the J2 satellite 
problem in the Delaunay-Similar elements, which is absolutely essential for fur- 
ther refinement in the Poincare'-Similar elements as presented here. 

The author also wishes to thank Stephen Starke of ACM, Inc. for checking 
the formulas and for generating the data in table I. 


This report was first issued as a Johnson Space Center Internal Note 
77-FM-52 JSC-13128 in November 1977. 


THE J2 PROBLEM IN DS4> VARIABLES 


This section briefly outlines the approximate solution (in the sense that 
only the short-period terms are eliminated) by the von Zeipel technique of 
the J2 problem in the DS4> variables. This outline is taken from reference 
4 where a more complete treatment of the J2 problem in the DS<j) variables may 
be found. 

The canonical DS<J> elements for the unperturbed case may be briefly 
described with the angle elements 


4> = true anomaly = t + constant 
g = argument of periapsis (classical w) 
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h = argument of the ascending node (classical 


and the action elements 
$ = associated with two-body energy 
G = total angular momentum 
H = axial component of the angular momentum 
L = the total energy 


The independent variable is T . 

The Hamiltonian for the oblateness perturbation J2 is 


is the perturbing Hamiltonian taken from reference 3. The following are abbrevi- 
ations for some functions of the DS(() elements: 


il = the time element 



y 

' _ T + constant 


F = F + eF, 
o 1 

F = - P/V2L 


( 1 ) 


o 


( 2 ) 


F^ is the unperturbed or two-body Hamiltonian in DS<1) elements and 



(1 + e cos <j>) 



(3) 



q = G - - $ + y/2V^ 

2 

2 2 2 

b = 1 - H /G = sin I (I is the inclination) 
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f = 


1 / 

pq 



( 4 ) 


The Hamiltonian equations of motion 


d4> 

3F 

d$ 

3F 

dT 


dT 

” H 

dg 

3F 

dG 

3F 

dT 

9G 

dT 

9g 

dh 

3F 

dH 

3F 

dT 


dT 

3h 

dH 

3F 

dL 

3F 

dT 


dT 

* dl 


are nonintegrable in this form. 

A transformation will now be found 




( 6 ) 


by means of the von Zeipel method. This method may be systematically employed 
to eliminate the short-period and long-period terms allowing the integration of 
the differential equations in the new primed variables to be obtained by simple 
quadrature. As mentioned above, for the present, only the short-period terms 
will be considered. 


The generating function S to be used to eliminating the short-period 
effects is 


S 


S + e S, 
o 1 


(7) 


where 
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= $'4> + G'g + H'h + L'^- 


( 8 ) 


is the identity transformation because 


9So 

$ = = $• 

d(J) 


9 So 

(|)' = — = (|) 
3$* 

3So 

G = - — = g' g' = = g 


H = 


L = 


3g 

9So 

3So 

32 , 


= H' 


= L' 


h' = 


3g* 

3So 

3H* 

3 So 
3L' 


= h 


= 2 , 


(9) 


S'] is also a function of the old (unprimed) angle elements and the new 
(primed) action elements. 


The function S] is a periodic function of the fast variable 4» and is 
found from the von Zeipel equations, 

Fo(f’,L') = Fo($’,L») (10) 


3Fq as] 

= -Fip($',G',H',L',(t),g) (11) 

3 $ ' 3 (|) 


Equation (10) shows that the unperturbed part of the Hamiltonian 
is left unchanged by the transformation. Equation (11) is used 
to determine the generating function Si. F-]p is the periodic portion of 
F-j and has a mean value of zero. The new perturbing Hamiltonian is found 
from 


f’ ($',G',H',L',_,g = F^(<K',G',H',L',<j),g) - F^p($' ,G' ,H' ,L» ,<j),g) (12) 

The periodic portion of F^ is found from equation (3) by inspection 
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e cos 4> 


=- 


1 

pq 


I ■ 

b Te 

~ 2 \_2 


cos (4) + 2g) + cos (2<|) + 2g) + - cos (3<t> + 2g) 


(13) 


Solving the partial differential equation (11) for S-j, 


Si 


= -(- f'e'b' + - f'e' sin 4 ) 
\2 , 3 / 


+ - f'e'b’ sin (4> + 2g) 


+ - f'b' sin (24> + 2g) 


+ — f'e'b' sin (34> + 2g) 
12 


(14) 


The abbreviations f, e', b', p', q' are the same as those in equa- 

tions (4) except they depend on the priined elements. 

The new (primed) Hamiltonian is 


F' = f' + £F' 
0 1 


(15) 


where from equations (3), (12), (13) 



( 16 ) 


Because F' depends only upon the elements G', H', L', the problem is 

now solved. 


The Hamilton equations of motion, equation (5) (but with primes on all the 
canonical variables) , yield the first-order solution. This operation will be 
performed in a following section. 
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The generating function (eq. (7)) where Sq and S-| are given by equa- 
tions (8) and (1M) defines the transformation between the primed and unprimed 
variables in either direction through first order. During this transformation, 
there are numerical problems for certain types of orbits. For both I = 

0 and e = 0, singularities occur because sin I and e occur in the denominators 
of certain terms. In the next section, the transformation to the PS<t> canonical 
elements will be given. The PS<}> elements do not have singularities for 
vanishing inclination and eccentricity. 

THE J2 PROBLEM IN THE PScfi ELEMENTS THROUGH FIRST ORDER 
The transformation among the DS4> and PScj) variables may be given as 

= (1) + g + h 

^ 2 ~ ~ - G) sin (g + h) 

= - Y2(G - H) sin h 
= 8. (17) 

A transformation between primed and unprimed PS(|) elements is sought. The 
approach taken will be to transform the generating function S-| given by equa- 
tion (14) from DS<1) elements to PS4> elements by direct substitution. That is, 

(17) 

Si($',G',H',L',4),g) ^Si(pJ,p^,p^,p^,a^,a2,d^, (18) 

Here it is assumed that equations (17) remain valid whether in the primed or 
unprimed system. 


p^ = "V2($ - G) cos (g + h) 
p^ = 'V2(G - H) cos h 
P4 = ^ 


The operation involved in making the transformation from S-| to S-| is 
done by temporarily dropping the distinction between the primed and unprimed PS4> 
variables. By direct substitution of equations (17) > a temporciry S-| (Pj^, dj^) is 
obtained in terms of the unprimed variables. Then, by making the assumption 
that the primed and unprimed PS4) variables are related through the order of 

the perturbing parameter e, that is Pi. = Pi!. + O(^) and a = a ' + 0(e), and 

eliminating Pj^ in favor of p^, the temporary Si becomes Si(Pj^ + 0(e), dj^) . 

Terms of 0(e) are then neglected to obtain the relation 


7 



The abbreviations of equation (4) may also be expressed in terms of the 
PS<t> variables by use of equation (17) 


P' 



+ P/V2P 


q’ = G' - 



+ 


y 



( 20 ) 


Also, there are now some new abbreviations: 



I 

s' = (G' + H')P^a^ 



( 21 ) 
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Numerically, through 0(e), it makes no difference in equations (19), (20), 
and (21) whether the primed or unprimed PS(j) variables are used. Algebrai- 

3S-| 3S-) 

cally in finding the derivatives — and , the difference between the 

3a, 3p * 

k 

primed and unprimed quantities is very important and must be strictly observed. 
For example, in the equation for G' (middle of equations (21) above), primed 
values of and p^ are used, but the unprimed value of is used. 


The unperturbed portion of the generating function is once again an iden- 
tity transformation. 


=0 "‘>2 °3 *'’■1 “4 


( 22 ) 


The generating function through first order is 


S + e 


(23) 


The transformation between the primed and unprimed variables is therefore 


3S 


3Si 


a, — — - — 0if + £ “ ' 

3p' 3p' 

k k 


3S 
k 3a, 


'k 


3Si 

3 ^ 


(24) 


The computation of the partial derivatives 


3Si 3S-1 

and 

9p’ ap, 

k ^ 


are straightforward and are facilitated by the introduction of additional 
abbreviations. Let 


9 


w = Q'/2p'q' 


3 



Y2 = sV2Q» 

Y3 = - ajc') 

= sin o ^ 

= sin 20^ 

= sin 30^ 

Cl = cos 0-| 

C2 = cos 20 1 

C3 = cos 30-1 (25) 

and the generating function (eq. (19)) becomes 

S-| = - — wy (26) 

q2 

It is at this point that the advantage of the PS<|) elements over the DS<t> 

3S^ 9S^ 

elements for the J_ problem occurs. In the derivatives r- and , 

2 ap^' 3p/ 


10 


no singularities occur for vanishing inclination and eccentricity. At the close 
of the second section, it was seen that singularities will occur for these two 
conditions in the DS<1) solution. However, when the Hamiltonian F' in the PS<t> 
equations is completed, it is found that the Hamiltonian differential equations 
are nonintegrable. However, the solution can be found by returning to the DS4> 
solution and transforming it by direct substitution into the PS<t> variables. 

This will be done in the following section. 


USE OF THE DS<1) SOLUTION 


From equations (5) (with primes on the DSt variables) and equation (16), 
the solution in the DS4> elements is given by 


(j)' = 
g' = 

h' = 


* i - 2/3) T - T 


- 2/3) . fb'] (t -tJ 

) 


h + - fb, T - T 

o 2 3 \ o 


Jl» = I' + 


o ,.3/2 


(2L«)- 


- %) 


+ - f^(b’ - 2/3 ) (t 


k - 


(27) 


where g^, h^, are the initial t = Tq values of the angle elements. 


The elements G', H*, L' are all constant and 


f = 


f: = 


f . = 


b, = 


b: = 


3f •/9'J*' 
9f /9G* 
9f /9L' 
9b */9 g' 
9b'/9H» 


( 28 ) 


1 1 


The approach to be taken here is to find expressions in terms of the DS4> 
elements that will relate the PS4> elements at an arbitrary value of the 
independent variable to the PS<|> elements at initiation. 


Add the first three equations of equation (27) to obtain 


<t>’ 


g» + h' = (})' + g* + h' + (t - T I 
® o o o \ Of 

+ ^ - 2/3) + f^Cb' - 2/3) + f'b^ + f'b^j (t - Tq) 


The left-hand side is dij by definition and the first term on the right is its 
initial value So 


a 


1 



+ ? [f’(b' - 2/3) + f^Cb' - 2/3) + f'b^ + f'b^ ] 



(29) 


Now add the second and third equations of equation (27) to obtain 


g* + h* = g; + h; + ^ V (t 



where 


V = f^^b' - 2/3)+ f'^b^ + b^) 


Taking the sine of both sides 


sin (g' + h') = sin ^g^ + h^ j cos || v^T - 

+ cos (g; + h;) sin v(t - 
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Now because = constant and G* = = constant, multiply both 

by - G' ) to obtain 


- G’ ) sin (g* + h') = \ - G^^ sin cc 

+ - °o) (so ^o) 


e / 

- vIt - 
2 \ 


II 


using equations (17) 


■ 


-V2$ 


” - G'^ sin ^ g' + h'j 



Finally 


cr* = o' cos 
2 20 


v/t - T \ 

- P„ sin 

5 v/'t - T \ 

[2 V 7 J 

20 

_2 V o/_ 


By a similar procedure 


» ^ » 


^2 = ^20 


a' = o'^ cos 


P3 = P30 


I i" - 



"r* / \“ 

1 + CJ* sin 

f v/t _ T I 

' 20 

2 V 0 ) 


I ^■'>3 - ^o) 

I - ^o) 


- P3O 


+ O^Q sin 


I ‘"“3 {; - ■'o) 

3 - \) 


2 3 


sides 



(30) 

(31) 

(32) 

(33) 
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f 

Because from equation (17) it follows from equation (27) that 


= °40 




A 3/2 

7 


(■' - ■'o) 5 fit' - 2^3] (t - ^o) O'!) 


And also from the fact that and L' are constant, from equation (17) 


p' =$'=$’= p 
o 


10 


p , = L’ = L' = p 


40 


(35) 

(36) 


The first-order DS<j) solution transformed by substitution into the PSij) solu- 
tion equation (29) through equation (36) provides the means of propagating the 

t f t ^ 

initial perturbed values, Pj^q and forward to their values Pj^ and 

at a specified (or otherwise determined) value of the independent variable. 

THE METHOD OF SOLUTION 


The unprimed and primed PS<|) elements are related through the generating 
function S(P',P) by equation (24). If, for example, the initial unperturbed 
elements have been obtained from, say, Cartesian coordinates, then the initial 
perturbed elements may be computed from 


3Si 

^ko ~ *^ko f 

9Si 

p' = p - e 

""ko ^ko aa,. 


k = 1 , ... ,4 


(37) 


Because the primed (perturbed) elements are related to the unprimed through 
order e, either the primed or the_unprimed values of Pj^ and Pj^ may be used 
in evaluating the derivatives of S-| on the right hand side of equation (37). 
The resulting error will be of order e^, which is beyond the accuracy of the 
first-order solution. 

Similarly, after a prescribed propagation interval, the unperturbed PS<|) 
elements may be computed again by using equation (24). 
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\C') = 



P^(T) 


P^(V * 



T 


k 


.,4 


(38) 


As in the initialization described in equation (37) either primed or unprimed 
values of and may be used in numerically evaluating the derivatives 

of the generating function in equation (38). 

The computation procedure can be visualized by referring to figure 1 . The 
conversion from Cartesian coordinates to PS<l5 elements and its inverse are 
given in the appendix. 


SPECIFIED TIME 


Because time is not the independent variable in the DS4> and PS<t> theories, 
provision must be made for propagating forward to a specified time. The proce- 
dure for accomplishing this is iterative. The following equations are evaluated 
with the perturbed (primed) values of and The prime is omitted for 

convenience. 

From the appendix, the time equation is 


t(T) = P^^ 




2 

e 


e sin 


) 


(39) 


A standard Newton iteration procedure is now used to solve for the value 
of T, which corresponds to the specified value of time, T 


T = 


^ 

(dt/dT)i 


(40) 


The derivative is given by 


dt 

dT 


/q 
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which is known at t.]. The iteration is continued until t - t-| is less than 
a previously specified tolerance. Each time the tolerance is not met, (the 

initial guess of t) is set equal to the last computed value of t from equation 
(40). 


The first value of the initial guess t-|, is chosen by referring to equation 
(27). To a good approximation 


T, - T z 

1 o 


T - 


l-(2l.')3/2 ^ I - 2/3) 


(41) 


The denominator of this equation is computed at the same time that the other 
quantities in the DS<t» solution are computed. 


SOME NUMERICAL RESULTS 


The solution for the J 2 problem in the singularity-free PS<}> variables 
has been programed on the Univac 1110 computer and is available for public use. 
This prototype program is called PSANS (for Poincare-Similar analytical short- 
period eliminated solution) . 


Several typical trajectories for satellites orbiting the Earth were 
computed for times corresponding to 1, 10, 50, and 100 revolutions of the satel- 
lites for each case. The initial conditions (in terms of classical orbital 
elements) are given, and the results of comparisons between the analytical solu- 
tion computed by this method (PSANS) and double-precision numerical integration 
are given in table I. It is seen that for a wide range of inclinations (0 to 
H/6 radians (0 to 30 degrees)) and eccentricities (0 to 0.1) the maximum errors 
are less than 1 kilometer. The errors were computed from the formula 


Error = 


- PSANS 


-INTEG 


at several points during each revolution. The maximum value of the error during 
a given revolution were output for the presentation in table I. 


CONCLUDING REMARKS 


From the numerical results presented, it can easily be seen that the 
analytical solution to the J 2 satellite problem in the PS(j) variables 
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overcomes the nuisance of the singularities due to small eccentricities and in- 
clinations that exist in the solution when Delaunay variables are used. 


The problem that remains is to systematically extend the solution to in- 
clude the cases ( 1 ) where both short- and long-period terms are included in the 
solution, (2) where the effects of the higher order zonal terms are included, 
and (3) where the effects of atmospheric drag are included. All of these new ex- 
tensions are currently being developed. The present version (PSANS), however, 
has demonstrated that it has a remarkable accuracy for medium lifetime satel- 
lites (less than 1 kilometer error in approximately 7 days). 
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TABLE I.- ANALYTICAL SOLUTION AND NUMERICALLY INTEGRATED SOLUTION 


COMPARING SIX TEST CASES 
(a) Initial conditions 


Orbital elements 

Case 

a 

e 

I 

(1) 

Q 

M 

1 

1 .050486 

0.015 

0 

TT/9 

n/9 

n/9 

2 

1 .050486 

0 

11/6 

n/9 

n/9 

n/9 

3 

1 .050486 

.001 

n/6 

11/9 

n/9 

n/9 

4 

1.160238 

.1 

11/6 

n/9 

n/9 

n/9 

5 

1 .070560 

0 

0 

0 

0 

0 

6 

1.070560 

0 

n/6 

0 

0 

0 


^The orbital elements are defined as folows: 

a = semimajor axis in Earth radii 

e = eccentricity 

I = inclination in radians 

w(or g) = argument of perigee in radians 

Jl!(or h) = argument of ascending node in radians 

M = mean anomaly in radians 
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TABLE I.- Continued 


(b) Revolution number, time function, maximum errors 


Rev 

Time , days 

Max error , m 

Case 1 

1 

0.06 

8 

10 

.63 

50 

50 

3.15 

254 

100 

6.30 

511 


Case 2 

1 

0.06 

11 

10 

.63 

100 

50 

3.15 

498 

100 

6.30 

995 


Case 3 

1 

0.06 

11 

10 

.63 

100 

50 

3.15 

498 

100 

6.30 

995 


Case 4 


1 

0.07 

10 

10 

.73 

84 

50 

3.66 

432 

100 

7.32 

872 


Case 5 

1 

0.06 

8 

10 

.65 

48 

50 

3.24 

239 

100 

6.48 

479 
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TABLE I.- Concluded 


(b) Concluded 


Rev 

Time, day 

Max error, m 

Case 6 

1 

0.06 

12 

10 

.65 

94 

50 

3.24 

470 

100 

6.48 

940 
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Input initial conditions 

and either t (final true 

anomaly) or T (final time). 
Convert initial conditions 

to p|^Q, (see appendix). 
Compute initial values of 

primed PS4> elements 

Pj^o) (eq.(24)). 

Compute DS<|) derivatives 

equations (27). 

Propagate to 

specified (or initial guess) 
^ (^s.j2^-(36))_^ 

Compute final values of 
unprimed PS(^ elements 
0 |^(t), P|^(t) equations (24). 


Compute trial value of t by 
Newton iteration. 


Update T. 

Convert PS 4 > elements to 
coordinates (see appendix). 



Figure 1.- Computation procedure for PSANS, 


APPENDIX 


COMPUTATIONAL PROCEDURE 

The equations that transform from position and velocity in Cartesian coordinates 
and the time to the PS<1> elements and the inverse transformation are listed 
below. 

1. From r, v, t, to p, a, (XTOPS). 

compute the perturbing potential V , 

G = £ X V (the angular momentum vector) 

G = |g| 



= -2G^/yj2(G + H) 
= -2G^/yj2(G + H) 
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q = -(p - Gj + - — L 

^ ^ 2 ^ 


p = - 


- G j + M/J2L 


Q = 


/r 


2 y / X 


1/2 




^2 = 


r-v 

r 


Q (v — + g] 


^2 = ^2 ^1 ~ 

P2 = Z2 sin cos 

e cos <t) = q(p2 cos sin 

e sin <j) = Q^p^ sin cos 

2 2 2 
e = (e cos (j>) + (e sin 4>) 


E - (j) = -2 tan 


-1 


e sin <|> 


= t 


(2L) 


3/2 


1 + Vi - + e cos <j) 

^E - 4> - - Vi - e sin <|) 
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2. From O, P to r, v, t (PSTOX) 



R = P- sin o + a cos o . 
j ' i 2 ' 



e cos 4> = Q ^p^ cos sin 

e sin <1) = Q ^p^ sin cos 

r = p/( 1 + e cos 4)) 

= r ^cos - a^R/2G^ 

X 2 = r ^sin - p^R/2G^ 

rR , 

X- = >(0 + H) 

^ zVg 

a = p, - 2 (p^2 , 


H = G - ! (pj 2 . <,32) 


• G / 
R = — (p 

r2 \ 


R = — ^P^ cos - 0^ sin 


in a,) 


^ (^2 ■*■ ^2 '^ l ) 


G ■ 


CsR 


X. = r ~ - rl — sin 0-] + 


r \^2 


2G 
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. ^2 

x_ = r — + r| — cos o-i - -i- 
? ” 2G 


= (rR + rR) 

2 2 2 
e = (e cos 4>) + (e sin 4>) 


E - (j) = -2 tan 


_1 


e sin 4> 


\1+Vl~6 + e cos (j) 


t = 0|j + 


(2P,) 


3/2 


E - ()) - 




e e sin (|) 
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